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Abstract 

The  overall  goal  of  this  research  was  to  develop  an  acoustic  wave  propagation 
model  using  well-understood  and  documented  computational  techniques  that  track  and 
quantify  an  air-borne  incident  acoustic  wave  propagated  around,  into  and  in  the  human 
head.  The  model  purpose  served  two  purposes:  (1)  determine  alternate  acoustic 
propagation  paths  to  the  cochlear  shell  that  exist  besides  the  normal  air-borne  acoustic 
propagation  path  (eardrum-ossical  path)  through  the  auditory  canal  and  (2)  quantify  sound 
pressure  amplitudes  in  the  cochlear  shell  (relative  to  the  air-borne  sound  pressure 
amplitude)  via  the  alternate  propagation  paths.  A  3D  finite-element  solid  mesh  was 
constructed  using  a  digital  image  database  of  an  adult  male  head.  Finite-element  analysis 
was  used  to  model  the  wave  propagation  through  the  fluid-solid-fluid  media.  Instantaneous 
acoustic  pressure  waveforms  were  recorded  at  various  positions  inside  and  outside  of  the 
head  model,  and  propagation  trajectories  (ray  paths)  were  constructed  and  evaluated  from 
wavefront  normals  as  a  function  of  frequency  and  incidence  angle.  The  acoustic  loss  across 
the  skull  was  estimated  to  be  *33  dB,  consistent  with  theoretical  estimates.  The 
computational  ray-path  results  and  the  theoretical  solutions  calculated  using  Snell’s  law 
gave  a  0.7°  difference  for  low-angle  incidence. 


Executive  summary  of  accomplishments 

1.  Validated  finite-element  analysis  (FEA)  general  processing  code  for  both  harmonic 
and  transient  solutions. 

2.  Constructed,  simplified  and  verified  transient  FEA  analyses  of  the  two-dimensional 
NIH  human  head. 

3.  Demonstrated  FEA  analysis  of  the  three-dimensional  Analyze  human  head. 

4.  Developed  a  ray  tracing  approach  to  graphically  and  quantitatively  represent  a 
three-dimensional  transient  process. 


Comprehensive  technical  report 
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report  titled  “Evaluation  of  Acoustic  Propagation  Paths  into  the  Human  Head”  and  authored 
by  W.  D.  O’Brien,  Jr.  and  Y.  Liu  is  attached.  In  addition,  Ms.  Y.  Liu  completed  her 
Electrical  Engineering  MS  thesis  in  2005  at  the  University  of  Illinois  at  Urbana-Champaign 
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on  this  research  topic.  Her  MS  thesis,  “Wave  Propagation  Study  Using  Finite  Element 
Analysis,”  as  well  as  appropriate  movie  files  (.avi),  are  available  for  downloading  at  the 
web  site:  www.brl.uiuc.edu/Downloads/. 
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Evaluation  of  Acoustic  Propagation  Paths  into  the  Human  Head 

Prof.  Dr.  William  D.  O’Brien,  Jr.  and  Ms.  Yuhui  Liu 

Bioacoustics  Research  Laboratory 
Department  of  Electrical  and  Computer  Engineering 
University  of  Illinois 
405  N.  Mathews 
Urbana,  IL  61801 
United  States  of  America 
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SUMMARY 

The  overall  goal  has  been  to  develop  an  acoustic  wave  propagation  model  using  well-understood  and 
documented  computational  techniques  that  track  and  quantify  an  air-borne  incident  acoustic  wave 
propagated  around,  into  and  in  the  human  head.  This  model  serves  as  a  computational  tool  to  elucidate  the 
acoustic  wave  propagation  around,  into  and  in  the  human  head.  Specifically,  the  model  determines  two 
features:  (1)  alternate  acoustic  propagation  paths  to  the  cochlear  shell  that  exist  besides  the  normal  air¬ 
borne  acoustic  propagation  path  (eardrum-ossical  path)  through  the  auditory  canal  and  (2)  sound  pressure 
amplitude  in  the  cochlear  shell  (relative  to  the  air-borne  sound  pressure  amplitude)  via  the  alternate 
propagation  paths.  A  3D  finite- element  solid  mesh  was  constructed  using  a  digital  image  database  of  an 
adult  male  head.  Coupled  acoustic-mechanical  finite-element  analysis  (FEA)  was  used  to  model  the  wave 
propagation  through  the  fluid-solid-fluid  media.  The  pressure  field  in  fluid  media  and  the  displacement  field 
in  solid  structures  were  computed  at  each  time  step.  Instantaneous  acoustic  pressure  waveforms  were 
recorded  at  various  positions  inside  and  outside  of  the  head  model,  and  propagation  trajectories  (ray  paths) 
were  constructed  and  evaluated  from  wavefront  normals  as  a  function  of frequency  and  incidence  angle.  The 
acoustic  loss  across  the  skull  was  estimated  to  be  approximately  33  dB,  consistent  with  theoretical  estimates. 
The  computational  ray-path  results  and  the  theoretical  solutions  calculated  using  Snell’s  law  gave  a  0.7° 
difference  for  low-angle  incidence;  10°  difference  was  obtained for  larger  angle  incidence. 


1.0  INTRODUCTION 

The  overall  objective  of  the  proposed  program  is  the  development  of  an  acoustic  wave  propagation  model 
using  well-understood  and  documented  computational  techniques  that  track  an  air-borne  incident  acoustic 
wave  propagated  around,  into  and  in  the  human  head.  This  model  will  serve  as  a  computational  tool  to 
elucidate  the  acoustic  wave  propagation  around,  into  and  in  the  human  head  to  evaluate  alternate  acoustic 
propagation  paths  that  exist  besides  the  normal  air-borne  acoustic  propagation  path  (eardrum-ossical  path) 
through  the  auditory  canal  to  the  cochlea,  a  cavity  in  the  temporal  bone.  Specifically,  the  computational 
model  will  evaluate  bone  conduction  that  is  produced  from  a  sound  field  in  air,  and  track  the  bone- 
conduction  acoustic  propagation  paths  into  the  cochlear  shell. 

The  validation  testing  of  the  computational  model  will  involve  human  subject  testing.  The  interplay  between 
human  subject  (experimental)  testing  and  theoretical  modeling  has  always  benefited  the  advancement  of 
science.  The  theoretical  model  provides  a  known  fundamental  basis,  and  the  experiment  provides  the 
evidence  to  evaluate  the  model.  When  the  experimental  and  modeling  results  essentially  agree,  then  a  basis 
is  provided  to  understand  at  a  fundamental  level  why  the  experiment  performs  in  the  way  it  does.  Further, 
and  most  important,  with  the  model  validated  against  experimental  observations,  the  model  then  becomes  a 
valuable  tool  to  test  other  hypotheses  without  the  need  to  conduct  expensive  and  time  consuming 
experiments. 
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2.0  MOTIVATION  AND  OBJECTIVE 

Noise-induced  hearing  loss  (NIHL)  has  been  an  important  issue  for  many  years.  Human  beings  that  are 
exposed  to  excessively  high  levels  of  noise  can  result  in  temporary  hearing  loss  or  permanent  hearing 
damage,  e.g.,  crews  that  must  work  near  military  aircraft  are  subject  to  severe  noise  environments,  with 
Sound  Pressure  Levels  (SPL)  reaching  145  to  150  dB.  Even  a  single  impulse  exposure  per  day  to  such 
intense  sounds  exceeds  permissible  noise  exposure  rules  [OSHA  1910.95],  Different  hearing  protection 
devices  (HPDs)  have  been  designed  to  prevent  NIHL,  such  as  earplugs,  earmuffs  and  helmets.  However, 
there  is  a  limit  to  how  much  attenuation  can  be  achieved  by  protecting  the  ear  canals  with  muffs  and/or 
plugs,  because  they  do  not  attenuate  sounds  that  are  conducted  to  the  inner  ear  by  the  hard  and  soft  tissues  of 
the  head  and  body. 

Bone-conduction  pathways  have  been  hypothesized  for  nearly  the  past  three-quarters  of  a  century  [Guild 
1936;  Barany  1938;  Wever  1954]  and  are  still  largely  unknown  [Sohmer  2000;  Sohmer  2004].  However, 
some  bone-conduction  pathways  have  been  proposed,  but  most  of  the  work  to  deduce  these  pathways  has 
been  via  direct  contact  methods  with  some  type  of  vibrating  device  [Wever  1954;  Freeman  2000;  Sohmer 
2000;  Stefan  2004],  Human  beings  sense  air-borne  sound  transmitted  by  bone  conduction  when  the  sound 
pressure  of  the  external  field  is  about  60  dB  above  the  threshold  required  for  hearing  via  the  eardrum-ossical 
path  [Bekesy  1948;  Zwislocki  1957].  This  suggests  that  noise  above  60  dB  SPL  cannot  be  fully  attenuated 
by  covering  or  blocking  the  auditory  canal. 

When  individuals  are  exposed  to  severe  noise  environments,  such  as  those  generated  by  aircraft  engines  and 
military  weapons  that  approach  and  even  exceed  1 50  dB  SPL,  even  if  hearing  protection  equipment  is  worn, 
they  may  be  subject  to  hearing  damage.  Furthermore,  NIHL  at  low  frequencies  (e.g.,  250  Hz  and  less)  are 
even  more  challenging.  Thus,  for  very  high  noise  environments  and  within  certain  frequency  ranges, 
conventional  HPDs  do  not  provide  sufficient  protection. 

Thus,  the  goal  of  the  proposed  program  is  to  deduce  bone-conduction  pathways  into  the  cochlear  shell.  From 
such  knowledge,  the  design  of  advanced  HPDs  can  then  be  based  on  well-founded  scientific  fundamentals;  it 
should  be  noted,  however,  that  the  program’s  goal  does  not  include  the  development  of  advanced  HPDs. 

The  normal  hearing  process  termed  “air  transmission”  or  “eardrum-ossical  path”  refers  to  air-borne  sound 
entering  the  auditory  canal  and  transduced  by  the  organ  of  Corti.  However,  other  propagation  path  processes 
are  considered  for  purposes  of  this  study.  The  process,  termed  “bone  conduction,”  refers  to  sound  entering  at 
any  location  other  than  the  auditory  canal  and  transduced  by  the  organ  of  Corti,  that  is,  other  than  the 
eardrum-ossical  path.  The  term  bone  conduction  does  not  necessary  imply  that  the  propagation  path  is 
entirely  bone.  The  distinction  between  these  two  processes  (eardrum-ossical  path  and  bone  conduction)  is 
necessary  because  the  acoustic  propagation  paths  are  evaluated  via  finite-element  analysis  (FEA)  procedures, 
and  the  eardrum-ossical  path  is  not  included  in  the  evaluation. 

For  purposes  of  this  study,  sound  levels  will  be  computationally  evaluated  to  and  in  the  region  referred  to  as 
the  cochlear  shell,  a  cavity  that  is  located  in  the  highly  porous  temporal  bone  that  contains  the  organ  of  Corti. 
The  cochlear  response  is  a  vectorial  integration  of  the  air-conduction  pathway,  different  bone-conduction 
pathways,  and  any  other  potential  alternative  pathways.  It  is  possible  that  each  of  these  pathways  is  more 
effective  at  different  frequencies.  Improving  the  current  hearing  protection  devices  thus  focuses  the  need  to 
better  understand  the  sound  propagation  around,  into  and  inside  of  the  human  head,  to  evaluate  different 
propagation  pathways  and  the  paths  that  are  taken  to  the  cochlea. 

3.0  FUNDAMENTALS  FOR  VALIDATION  STUDIES 

3.1  Wave  Equation 

The  necessary  assumptions  used  in  ANSYS  acoustic  analysis  include: 


3 


•  The  fluid  is  compressible,  but  only  relatively  small  pressure  changes  with  respect  to  the  mean 
pressure  are  allowed. 

•  The  fluid  is  inviscid  (no  viscous  dissipation). 

#  There  is  no  mean  flow  of  the  fluid. 

#  The  mean  density  and  mean  equilibrium  pressure  are  uniform  throughout  the  fluid. 


The  acoustic  wave  equation  is  given  by  [Pierce  1989] 


c^P 

at2 


c  V  P  =  0 , 


(1) 


where  P  is  the  acoustic  pressure,  c  is  the  sound  speed  in  the  fluid  medium,  t  is  time  and  V  is  the  Laplacian 
operator. 

The  finite-element  statement  of  the  wave  equation  is  given  by 

tMe]{Pe}  +  [KeP]{Pe}  +  Po[Re]TK}  =  {0},  (2) 

where 

[Me  ]  =  -j  J  d(vol)  is  the  fluid  mass  matrix  (2a) 

[Kep]  =  J  [B]T[B]d(vol)  is  the  fluid  stiffness  matrix  (2b) 

Po  [Rj  =  Po  J {N}{n}T  {N  '}t  d(S)  is  the  fluid-structure  coupling  mass  matrix  (2c) 

s 

{N}  is  the  element  shape  function  for  pressure 
{N '}  is  the  element  shape  function  for  displacements 
{Pc}  is  the  nodal  pressure  vector 

{ue}  =  {uxe},  {uye},  {uze}  is  the  nodal  displacement  component  vectors 
{n}  is  the  unit  normal  at  the  fluid  boundary 
vol  is  the  volume  of  domain 

S  is  the  surface  where  the  derivative  of  pressure  normal  to  the  surface  is  applied 

Well-characterized  analytical  solutions  [Morse  1968]  are  compared  with  the  finite-element  simulation  results 
to  validate  the  analysis. 

3.2  Scattering  by  Rigid  Cylinder  (2D  Rigid  Cylinder  Case) 

The  acoustic  scattering  from  a  rigid  cylinder  of  radius  a  due  to  a  plane  wave  propagating  in  a  direction 
perpendicular  to  the  cylinder’s  axis  is  given  by 


pinc  =  P0eik(rcos*-ct)  =  P0[J0(kr)  +  2]T  im  cos(m(|>)Jn(kr)]e- 

m=l 

Psca  =  X  Am  cos(m<|))[Jm(kr)  +  iNin(kr)]e_i“t 


m-0 


Am=-emP0im+1e-'^  sinYm 
tanYm  = 

°  Nj(ka)  Nm+1(ka)-Nm_1(ka) 


(3) 
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where 


pinc  is  the  incident  plane  wave 
psca  is  the  scattered  wave 
a  is  the  cylinder  radius 

r  is  the  radial  distance  from  the  center  of  the  cylinder 

<|>  is  the  axial  angle  for  cylindrical  coordinates 

c  is  the  speed  of  sound  in  the  medium  surrounding  the  cylinder 

k  is  the  oo/c  =  2k/X,  wavenumber 

ym  is  the  phase  shift 

em  =  1  when  m  =  0;  £m  =  2  when  m  >  0 

The  total  acoustic  pressure  at  the  surface  of  the  cylinder  at  an  angle  <|>  relative  to  the  x  axis  is  given  by 


4P, 


V  COS(m(|)) 


Pa  =  Pine  +  Psca  = 

rcka  m=0  E 


where  Em,  the  radiation  amplitude  for  a  cylinder,  is 
E  E 


8 

!27tka 

when  ka  »  m  + 1  /  2 ,  and 
4 


m>0 


2 

Ttka 


Jtka 


m!  2  , 

’  *im>0  ~  2tc 'ka  ’ 


(4) 


(4a) 


(4b) 


when  ka  «  m  + 1  /  2 . 

3.3  Scattering  by  Rigid  Sphere  (3D  Rigid  Sphere  Case) 

The  acoustic  scattering  from  a  rigid  sphere  of  radius  a  due  to  a  plane  wave  propagating  to  the  right  along  the 
polar  axis  is  given  by 


Pi„c  =  P0eik(rcos»-ct)  =  P0  £  (2m  +  l)imPm  (cos,d)jm  (kr  )e-tot 

m=0 

Psca  =  -PoI(2m  +  1)re-i8m  sin 5mPm  (cos  $)[  jm  (kr)  +  inm  (kr)]e'i<0‘  (5) 

m=0 

where  psca  is  the  wave  scattered  from  the  sphere  whose  center  is  the  polar  origin.  The  total  acoustic  pressure 
at  a  point  on  the  sphere  an  angle  $  relative  to  the  polar  axis  (note  that  the  point  0  =  0°  is  the  location 
farthest  away  from  the  source  of  the  sound)  is  given  by 


Pa  =  (Pine  +  Psca)r=a  =  Po(ka)'2  £ 


m=0 


(2m +  1) 

B„ 


Pm(cosfl)e 


-i(Sm  -trnn-cot) 


~  (1  + 1  ika  cos  'd)P0e_,rat ,  ka  « 1 


(6) 
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3.4  Scattering  by  Nonrigid  Sphere  (3D  Elastic  Sphere  Case) 

The  scattering  from  a  sphere  that  is  not  perfectly  rigid  is  more  complicated.  The  total  acoustic  pressure  at  a 
point  on  the  sphere  an  angle  $  relative  to  the  polar  axis  (note  that  the  point  A  =  0°  is  the  point  farthest 
away  from  the  source  of  the  sound)  is  given  by 


Pa  =  Po X (2m  +  l)imPm(cosd)[jm(kr) - i(l  +  Rm)hm (kr)] , 

m=0 


where 


Rm,  the  reflection  coefficient,  is  1  +  Rm  =  2 


(3m,  the  effective  admittance,  is  (3ra  =  i 


;  PC 


PeCe 


jm(ka)  +  iftmjm(ka) 
hm(ka)  +  i(5mhm(ka) 

Jm(M) 

Jm(kea) 


ke  =  co/ce 

p,  c  are  the  density  and  speed  of  the  medium  surrounding  the  sphere 
pe,  ce  are  density  and  speed  of  the  sphere  material 


4.0  VALIDATION  STUDIES 


(7) 


(7a) 

(7b) 


The  ANSYS  (ANSYS,  Inc.,  Canonsburg,  PA)  finite-element  analysis  (FEA)  general  processing  code  for  both 
harmonic  and  transient  solutions  has  been  validated  using  well-understood  two-dimensional  (2D)  and  three- 
dimensional  (3D)  models.  Harmonic  (continuous-wave)  validation  studies  were  conducted  for  (1)  2D  rigid 
cylinder,  (2)  3D  rigid  sphere  and  (3)  3D  elastic  sphere.  Transient  (one-cycle  sine  wave)  validation  studies 
were  conducted  for  (1)  homogeneous  air,  (2)  2D  elastic  cylinder  shell,  and  (3)  3D  elastic  (water)  sphere.  For 
all  of  the  validation  studies,  either  a  harmonic  or  transient  acoustic  plane  wave  was  initiated  in  air.  The  air¬ 
borne  acoustic  wave  was  incident  on  the  geometric  model  that  was  either  rigid,  elastic  or  water.  Water  is  an 
ideal  fluid  for  this  study  because  it  has  acoustic  propagation  properties  similar  to  those  of  brain  and  other  soft 
tissues.  In  all  cases,  the  computational  solutions  of  acoustic  pressure  distribution  agreed  well  with  the 
analytic  solutions  of  acoustic  pressure  distribution. 

To  be  familiar  with  different  types  of  acoustic  analysis  in  ANSYS  and  ensure  building  our  FEA  model  in  a 
proper  way,  several  case  studies  were  carried  out  on  different  geometric  models.  For  these  validation 
studies,  the  simulation  results  are  compared  with  the  theoretical  solutions  to  evaluate  the  model.  Herewith 
are  the  symbols  used: 

f  =  center  frequency  of  the  incident  wave  (kHz) 
c  =  sound  speed  in  the  medium  (m/s) 

X.  =  wavelength  in  the  medium  (m) 
p  =  material  density  (kg/m3) 

a  =  cylinder  radius  /  shell  outer  radius  /sphere  radius  (3D)  (m) 
d  =  shell  thickness  (m) 

r  =  radial  distance  from  the  center  of  the  cylinder 
BOUND  =  FEA  absorbing  boundary  radius  (m) 

(pinc  =  incident  angle,  the  angle  between  incident  wave  and  +x  axis  (degree) 

Xinc  =  incident  wave  position  (m) 

Pinc=  incident  wave  pressure  amplitude  (Pa) 
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T  =  period  of  the  incident  wave 

DPW  =  number  of  mesh  divisions  per  wavelength 

ITS  =  integration  time  step  size 

4.1  2D  Rigid  Cylinder  (Harmonic  Analysis  Case) 

Figure  1  shows  the  FEA  model  used  in  ANSYS  for  the  2D  rigid  cylinder  case  (§3.2).  The  air-borne  incident 
harmonic  plane  wave  propagates  in  the  +x-axis  direction.  When  the  propagating  wave  meets  the  rigid 
cylinder  target,  acoustic  pressure  is  distributed  around  the  cylinder  surface. 


Parameters: 
f  =  3  kHz 
cair  =  340  m/s 

*air  =  340/ 3000  =  0.1 13  m 
pair=  1.2  kg/m3 
a  =  0.4  Xair  =  0.0452  m 
BOUND  =  1.4  Xair  =  0.158  m 

tpinc  0 

Xinc=  -0.625  Aair  =  0.0706  m 
Pine  =  1  Pa 

Fig  1 :  The  2D  rigid  cylinder  FEA  model  DPW  =  20 

generated  by  ANSYS.  Red:  air;  Light 

blue:  rigid  cylinder;  Dark  blue:  connecting 

elements;  Yellow:  location  of  plane  wave 

initiation. 

Figure  2  shows  the  ANSYS  simulation  results  (circles)  of  the  pressure  distribution  on  the  cylinder  surface 
compared  with  the  analytical  solutions  (line),  and  along  the  radial  component.  Good  agreement  is  obtained. 


total  pressure  on  cylinder  surface,  n  »  0.4* 


total  pressure  along  ratftus.  <> «  0.  a  ■  <1.4* 


fl\ 


Fig  2:  2D  rigid  cylinder  simulation  results  (circles)  vs.  analytical  solutions  (lines).  Top  plot:  total  pressure  on  cylinder 
surface  (note  that  the  point  <j)  =  0”  is  the  point  farthest  away  from  the  source  of  the  sound).  Bottom  plot:  total  pressure 
along  radius,  observation  angle  <J)  =  0°  (note  that  r  =  radial  distance  from  the  center  of  the  cylinder). 
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4.2  3D  Rigid  Sphere  (Harmonic  Analysis  Case) 

Figure  3  shows  the  FEM  model  used  in  ANSYS  for  the  3D  rigid  sphere  case  (§3.3).  The  air-borne  incident 
harmonic  plane  wave  propagates  in  the  +x-axis  direction.  This  case  is  of  considerable  practical  importance 
because  many  scattering  objects  are  more  or  less  spherical.  When  the  propagating  wave  meets  the  rigid 
spherical  target,  acoustic  pressure  is  distributed  around  the  spherical  surface. 


Parameters: 
f  =  3  kHz  ' 

Cair  =  340  m/s 

Xair=  340/ 3000  =  0.1 13  m 
pair=  1.2  kg/m3 
a  =  0.4  lair  =  45.2  mm 
BOUND  =  2.4  1^  =  0.271  m 

tpinc  0 

Xinc=  -2.0  lair  =  0.226  m 
Pine  =  1  Pa 
DPW  =  10 


Fig  3:  The  3D  rigid  and  nonrigid  sphere  FEA  model. 

Figure  4  shows  the  ANSYS  simulation  results  (circles)  of  the  pressure  distribution  on  the  spherical  surface 
compared  with  the  analytical  solutions  (line),  and  along  the  radial  component.  Good  agreement  is  obtained. 
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Fig  4:  3D  rigid  sphere  simulation  results  (circles)  vs.  analytical  solutions  (lines).  Top  plot:  total  pressure  on  the  surface 
of  the  center  sphere  cross  section  (note  that  the  point  $  =  0'  is  the  point  farthest  away  from  the  source  of  the  sound). 
Bottom  plot:  total  pressure  along  +z  axis,  observation  angle  <j)  =  0°  (note  that  r  =  radial  distance  from  the  center  of  the 
sphere). 


4.3  3D  Elastic  Sphere  (Harmonic  Analysis  Case) 


For  the  perfectly  rigid  cases  (§4.1  and  §4.2)  the  propagated  wave  does  not  enter  the  object.  However,  for 
real  objects,  the  propagated  acoustic  wave  will  enter  the  object.  Figure  3  shows  the  FEA  model  used  in 
ANSYS  for  the  3D  nonrigid  sphere  case  (§3.4).  Here  the  nonrigid  sphere’s  density  and  speed  are  assumed  to 
be  twice  those  of  air  (parameters  same  as  in  Fig  3  except  for  csphere  =  680  m/s  and  psphere  =  2.4  kg/m3).  The 
air-borne  incident  harmonic  plane  wave  propagates  in  the  +x-axis  direction.  This  case  is  of  considerable 
practical  importance  because  many  scattering  objects  are  more  or  less  spherical.  When  the  propagating  wave 
meets  the  nonrigid  spherical  target,  acoustic  pressure  is  distributed  around  the  spherical  surface,  and  also 
enters  into  the  sphere. 

Figure  5  shows  the  ANSYS  simulation  results  (circles)  of  the  pressure  distribution  on  the  elastic  spherical 
surface  compared  with  the  analytical  solutions  (line),  and  along  the  radial  component.  Fair  agreement  is 
obtained. 
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1  - - - T — — - . r- - T — — — • - t - - - “T 


©  ANSYS 
Theoretical 


0  0.2  0.4  0-6  0.8  1  1.2  1-4 

r/X 


Fig  5:  3D  elastic  sphere  simulation  results  (circles)  vs.  analytical  solutions  (lines).  Top  plot:  total  pressure  on  the 
surface  of  the  center  sphere  cross  section  (note  that  the  point  <j>  =  0°  is  the  point  farthest  away  from  the  source  of  the 
sound).  Bottom  plot:  total  pressure  along  +z  axis,  observation  angle  point  <j>  =  0°  (note  that  r  =  radial  distance  from  the 
center  of  the  sphere). 


4.4  2D  Analysis  in  Homogenous  Media  (Transient  Analysis  Case) 

Transient  analysis  is  a  more  realistic  approach  than  harmonic  analysis  for  modeling  the  incident  acoustic 
pressure  waveforms  onto  and  in  the  human  head.  However,  transient  analysis  is  more  involved  because  it 
requires  added  computer  resources  and  more  of  our  resources,  in  terms  of  the  “engineering”  time  involved. 
To  save  a  significant  amount  of  these  resources,  preliminary  studies  were  conducted  to  understand  further  the 
physics  of  the  problem  for  validation  purposes,  that  is,  analyzing  a  simpler  model  provides  better  insight  into 
the  problem  at  minimal  cost. 


9 


The  2D  transient  analysis  (air  is  treated  as  the  object)  is  evaluated  to  determine  whether  the  computation 
domain  space  introduces  artifacts.  Here,  a  pulse  wave  (gated  sinusoid  wave)  initiated  in  homogenous  air 
medium  is  simulated.  A  simple  geometry  model  (Fig  6)  is  used  that  is  similar  to  the  2D  cylinder  harmonic- 
analysis  model  except  that  the  2D  cylinder  is  air.  The  purpose  of  using  an  air-filled  cylinder  submerged  in 
air  is  to  consider  possible  artifacts  that  might  be  caused  by  different  meshes.  The  pulse  is  applied  in  air  and 
propagates  in  the  +x  direction.  The  locations  A  (-0.9A,  0),  B  (-0.65a,  0),  C  (0.4A,  0),  H  (0,  0.4/1.)  (in  Fig  6) 
denote  where  the  acoustic  pressure  waveforms  are  plotted  in  Fig  7  to  evaluate  mainly  DPW  and  ITS. 

Parameters: 
f  =  3  kHz 

Air  (Cair  =  340  m/s,  Aajr =  0. 1 1 3  m,  =  1.2  kg/m3) 

(Pine  =  0° 

Pine  =  1  Pa 
DPW  =  variable 
ITS  =  variable 
a  =  0.4  =  45.2  mm 

BOUND  =  a  +  0.9  Aajr  or  a  +  1 .4 


Fig  6:  The  wave  propagation  through  air  only. 
Locations  A,  B,  C  and  H  are  denoted  as  red  letters 
where  A,  B  and  C  are  on  the  x  axis  and  H  is  on  the 
y  axis. 


Fig  7:  Pressure  waveforms  at  four  locations  in  homogeneous  air  medium  to  test  DPW  and  ITS. 

DPW  (number  of  mesh  divisions  per  wavelength,  affects  spatial  resolution)  and  ITS  (integration  time  step 
size,  affects  temporal  resolution)  are  two  key  transient-analysis  parameters,  and  BOUND  (FEA  absorbing 
boundary)  determines  the  degree  to  which  the  computational  domain  can  be  assumed  infinite  in  extent. 
Note:  the  incidence  plane  is  located  left  of  the  center  plane  in  the  computation  domain  and  BOUND  affects 
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the  degree  of  the  incident  wave’s  distortion  along  the  circular  boundary  of  the  computation  domain.  In 
general,  the  mesh  must  be  fine  enough  to  resolve  the  largest  dominant  frequency.  A  general  guideline  is  to 
have  at  least  20  elements  per  wavelength  along  the  propagation  direction,  that  is,  DPW  >  20  [ANSYS  6.1 
Documentation:  Modeling  and  Meshing  Guide],  For  transient  analysis,  using  at  least  twenty  points  per  cycle 
at  the  center  frequency  yields  a  reasonably  accurate  solution.  ITS  is  given  by  ITS  <  T/DPW,  where  T  =  1/f,  f 
is  the  center  frequency  of  the  incident  pulse  wave  [ANSYS  6.1  Documentation:  Transient  Dynamic 
Analysis],  The  acoustic  absorbing  boundary  should  be  located  at  least  0.2A  from  the  target  object  [ANSYS 
6.1  Documentation:  Coupled-Field  Analysis  Guide],  Six  combinations  of  DPW,  ITS,  and  BOUND  were 
evaluated  (Fig  7).  It  was  observed  that,  in  general,  larger  BOUND  enhanced  the  computational  accuracy 
(Tests  5  and  6  vs.  Tests  1,  2,  3,  and  4)  while  larger  DPW  and  smaller  ITS  yield  insignificant  improvements 
(Test  2  vs.  Test  1,  Test  3  vs.  Test  2).  For  future  transient-analysis  cases,  DPW  =  20  and  ITS  =  TOO  will  be 
considered  to  be  sufficient  and  BOUND  will  be  further  evaluated  in  next  section. 

4.5  2D  Elastic  Cylinder  Shell  (Transient  Analysis  Case) 

Table  1  lists  the  properties  for  most  components  of  human  head  [Goss  1978;  Goss  1980;  Duck  1990]. 
Among  these  properties,  sound  speed  in  air  is  the  lowest  and  thus  the  wavelength  in  air  is  the  shortest.  Thus, 
for  the  FEA  model,  the  smallest  elements  used  are  based  on  propagation  in  air.  In  the  other  words,  to 
achieve  good  computation  resolution  for  the  mixed-property  model,  the  optimal  parameters  for  simulating 
propagation  in  air  is  determined. 


Table  1 :  Material  Properties  in  the  Human  Head 


Material 

Speed  of  Sound 
(m/s) 

Densities 

(kg/m3) 

Air 

340 

1.2 

Water 

1500 

1000 

Soft  Tissues 

1520-1580 

980-1010 

Lipid-based  tissues 

1400-1490 

920-940 

Collagen-based  tissues 

1600-1700 

1020-1100 

Aqueous  humor 

1002-1006 

1500 

Vitreous  humor 

1090 

1530 

Blood 

1580 

1040-1090 

Brain  -  grey 

1532-1550 

1039 

Brain  -  white 

1043 

Skull  -  compact  inner  and  outer  tables 

2600-3100 

1900 

Skull  -  Spongy  diploe 

2200-2500 

1000 

A  transient  analysis  on  a  2D  soft  (elastic)  cylinder  shell  is  evaluated  that  includes  property  values  of  typical 
tissue  (Table  1).  Figure  8  shows  the  2D  soft  shell  geometry  used  in  the  simulations.  In  this  case  study,  the 
effect  of  BOUND  is  further  evaluated  in  object  spaces  of  different  sizes. 

Two  2D  cylinder  shell  tests  were  carried  out  with  different  FEA  absorbing  boundary  on  the  same  target 
object  (Table  2). 


Table  2:  Small  2D  Shell  Transient  Analysis 


Parameters 

Test  1 

Test  2 

Shell  inner  radius  (a) 

0.1  h* 

0.1  A-air 

Shell  thickness  (d) 

mmmmm 

BOUND 

Incident  wave  position 

-  (a  +  d  +  0.2)  Aajr 

-  (a  +  d  +  1.2)  Aair 

Distance  between  B  and  E  (Fig  8) 

0.57  Aair 

2.56  Aair 
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Parameters: 
f  =  3  kHz 

Outer  medium:  air  (cair  =  340  m/s, 

Kk  =  o.  1 1 3  m,  pair  =  1 .2  kg/m3) 

Inner  medium:  water  (cwater=  1500  m/s, 
Pwater  =  1000  kg/m3) 

Soft  Shell:  (cSheii  =  2800  m/s, 
p shell  =  1900  kg/m3) 
tPinc  —  H 
Pj„c  =  1  Pa 
DPW  =  30 
ITS  =  T/30 

Computational  length  =  5  T 
BOUND  =  variable 

Fig  8:  The  soft  shell  2D  FEA  model  generated  in  ANSYS. 


Figure  9  shows  x-axis  acoustic  pressure  waveforms  at  locations  A,  B,  C  and  D  (Fig  8).  The  circular  shape  of 
the  absorbing  boundary  causes  distortion  of  the  propagated  wave  along  the  boundary,  and  thus  induces 
additional  artificial  incident  waves  along  non  +x-axis  direction.  When  the  absorbing  boundary  is  positioned 
at  a  greater  distance  from  the  object,  less  undesirable  interference  reaches  the  target  object  in  the  same 
computational  length  scale.  The  distance  between  B  and  E  (E  is  at  the  absorbing  boundary;  Fig  8)  in  Table  2 
provides  a  rough  estimate  of  the  shortest  possible  time  required  for  the  interference  wave  to  reach  B  from  the 
absorbing  boundary.  Comparing  the  first  period  of  the  pressure  waveform  at  B,  the  waveform  in  Test  2  is 
more  symmetric  than  in  Test  1.  A  reasonable  explanation  is  that  it  requires  about  0.5T  for  the  interference 
wave  at  E  to  reach  B  in  Test  1  whereas  it  is  requires  about  2.5T  in  Test  2.  It  is  also  suspected  that  the 
significant  variation  of  the  pressure  waveform  shape  at  C  in  Test  1  at  about  IT  after  the  incident  wave 
reaches  C  is  due  to  the  interference  wave  from  the  absorbing  boundary  reaching  the  target  shell  cylinder  by 
that  time.  In  Test  2,  this  interference  is  smaller,  but  still  present,  because  it  takes  longer  for  the  interfering 
waves  at  a  further  absorbing  boundary  location  to  reach  the  shell  cylinder. 


Fig  9:  Acoustic  pressure  waveforms  in  Test  1  (a,  left)  and  Test  2  (b,  right) 

To  further  demonstrate  this  observation,  two  large  2D  shell  cylinder  tests  were  carried  out  with  different  FEA 
absorbing  boundary  (Table  3).  In  Test  3,  the  absorbing  boundary  is  located  the  same  distance  from  the  shell 
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cylinder  as  in  Test  1.  And  in  Test  4,  the  absorbing  boundary  is  moved  further  away  than  in  Test  3  but  not  as 
far  as  in  Test  2. 


Table  3:  Large  2D  Shel 

Transient  Analysis 

Parameters 

Test  3 

Test  4 

Shell  inner  radius  (a) 

0.58  Aair 

0.58  A^air 

Shell  thickness  (d) 

0.09  Lajr 

0.09  Aair 

BOUND 

Incident  wave  position 

HEQEEEfB 

-  (a  +  d  +  0.4) 

Distance  between  B  and  E  (Fig  8) 

0.81  A,air 

1.34  A,air 

Figure  10  shows  x-axis  acoustic  pressure  waveforms  at  locations  A,  B,  C  and  D  (Fig  8).  In  Test  4,  an 
abnormal  peak  pressure  amplitude  occurs  at  C  about  IT  after  the  incident  wave  reaches  this  position,  which 
is  not  observed  in  the  other  tests.  Again  it  is  very  likely  due  to  the  interfering  waves  coming  from  the 
absorbing  boundary.  Similar  to  the  previous  small  shell  cylinder  case,  more  symmetric  pressure  waveforms 
for  the  first  period  are  found  in  Test  4.  In  both  Tests  3  and  4,  significant  variation  of  the  pressure  waveform 
shape  at  B,  C  and  D  are  observed  about  IT  after  the  incident  wave  reaches  that  position,  which  implies  the 
absorbing  boundary  is  not  sufficiently  further  away  from  the  shell  cylinder.  Based  on  the  observations  (§4.2 
and  §4.3),  it  is  concluded  that  FEA  absorbing  boundary  has  a  significant  impact  on  the  computation  accuracy 
and  larger  BOUND  is  preferred  to  decrease  the  interference  due  to  the  circular  absorbing  boundary  and  to 
better  simulate  the  incident  plane  wave. 


Fig  10:  Acoustic  pressure  waveforms  in  Test  3  (a,  left)  and  Test  4  (b,  right) 

4.6  3D  Elastic  (Water)  Sphere  (Transient  Analysis  Case) 

In  this  case  study,  transient  analysis  is  conducted  on  a  3D  sphere  of  water  that  is  submerged  in  air  (§3.4),  as 
shown  in  Fig  11.  The  acoustic  pressure  waveforms  are  shown  at  different  locations  (Fig  1 1)  in  Fig  12.  The 
locations  are  A  (0,  0,  -0.9A),  B  (0,  0,  -0.4A.),  C  (0,  0,  0),  D  (0,  0,  0.4 A,).  From  the  pressure  waveforms  it  is 
determined  that  about  0.5  T  is  required  for  the  incident  wave  to  travel  from  A  to  B  whereas  the  travel  time  is 
around  0. 1  T  from  B  to  C.  These  travel  time  estimates  agree  with  theory  and  thus  validates  the  temporal 
accuracy  of  the  FEA. 
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Parameters: 
f  =  3  kHz 

Air:  cair  =  340  m/s,  X,air  =  0.1 13  m,  pair  =  1.2  kg/m3 
Water  sphere:  Cwater  =  1500  m/s,  pwater  =  1000  kg/m' 
a  =  0.4  Xair 
BOUND  =  1.3  X,* 

Xi„c  =  0.9  Kr 

tPinc  —  0 

Pine  =  1  Pa 
DPW  =  20 
ITS  =  T/20 


Fig  1 1 :  The  3D  water  sphere  FEA  model. 


Furthermore,  the  instantaneous  acoustic  intensity  was  calculated  based  on  the  FEA  acoustic  pressure  data 
from 


PA(t)2 

Pair^air 


IC  = 


Pb(02 

)  C 

'water  water 


(8) 


where  p(t)  are  the  peak  values  of  the  acoustic  pressure  at  their  respective  locations.  The  instantaneous 
acoustic  intensity  insertion  loss  (IL)  from  air  to  water  was  estimated  from 


ILdB  ~  101og10 


(\  \ 

Apeak,C 

^  ^peak,A  ) 


(9) 


From  the  FEA  data  (Fig  12),  the  intensity  loss  from  air  to  water  was  estimated  to  be  -34  dB,  a  reasonable 
value  relative  to  a  theoretical  loss  estimate  (-33  dB)  for  normal  incidence  at  a  planar  air- water  boundary. 


f=3kHz 


Fig  12:  Acoustic  pressure  waveforms  for 
the  water  sphere  case  at  locations  A,  B, 
C,  and  D  (Fig  11). 
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5.0  HUMAN-HEAD  RELATED  STUDIES 
5.1  NIH  Human  Head 

Transient  FEA  procedures  of  the  2D  NIH  human  head  (National  Library  of  Medicine’s  Visible  Human 
Project,  National  Institutes  of  Health,  Bethesda,  MD)  have  been  constructed,  simplified  and  verified.  The 
male  anatomic  dataset  was  used.  The  anatomic  dataset  consists  of  0.33-mm-wide  transverse  sections  of  the 
head,  with  each  section  2048  pixels  by  1216  pixels  and  each  pixel  8-bit  RGB  scale.  A  simplified  2D  human 
head  analysis  was  conducted  using  one  of  the  0.33-mm-wide  transverse  sections.  The  first  challenge 
(computational,  not  scientific)  was  to  convert  the  NIH  human  head  digital  raster  image  dataset  into  a  vector 
format  file  for  import  to  ANSYS  that  uses  IGES-format  vector  data.  Then,  the  outer  surface  of  the  head  was 
segmented  to  yield  only  the  head  contour.  Skull  was  modeled  as  a  1-cm-thick  layer  immediately  inside  the 
human  head  contour.  The  interior  region  was  modeled  as  water  (Fig  13). 

Transient  (one-cycle  sine  wave)  analyses  were  performed  at  4  frequencies  (0.125,  1,  3  and  10  kHz)  and  two 
incident  angles  (0°:  towards  the  right  side,  and  45°:  approximately  towards  the  right  cheek)  (Fig  14). 


Fig  13:  Development  of  two-dimensional  geometric  model  of  human  head  based  on  anatomic  image. 

The  instantaneous  acoustic  pressure  waveforms  were  recorded  at  4  locations,  A,  C,  F  and  H  (Fig  15),  which 
are  at  approximately  the  left  and  right  ear  locations  near  each  side  of  the  skull,  and  shown  in  Fig  16. 
Validation  of  the  process  was  accomplished  by  estimating  the  acoustic  intensity  loss  across  the  skull  from  the 
instantaneous  acoustic  pressure  waveforms  on  each  side  of  the  skull  (one  location  in  air  and  the  other 
location  in  water).  The  acoustic  loss  (Eqns  8  &  9;  §4.5)  across  the  skull  was  estimated  from  the  FEA  data  to 
be  approximately  26  dB,  reasonably  consistent  with  theoretical  estimates  (33  dB)  considering  this  is  a  2D 
analysis.  The  instantaneous  acoustic  intensities  are  shown  in  Fig  17  for  two  tests  with  same  incident 
frequency  (3  kHz)  and  the  two  incident  angles  (0°  and  45°). 

5.2  Analyze  Human  Head 

The  transient  FEA  procedure  of  the  3D  Analyze  (Mayo  Clinic,  Rochester,  MN)  human  head  model  has  also 
been  demonstrated.  The  human  head  model  was  derived  from  the  3D  Analyze  MRI  dataset.  The  Analyze 
human  head  outer  surface  was  segmented  slice-by-slice  as  shown  in  Fig  18. 
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Parameters: 

Outer  medium:  air  (cajr  =  340  m/s,  X^  =  0.1 13  m, 
pair=  1.2  kg/m3) 

Inner  medium:  water  (cwater =  1500  m/s,  pwater  =  1000  kg/m3) 
Human  skull  (p  =  1412  kg/m3,  Young’s  modulus  =  6.5  GPa, 
Poisson’s  ratio  =  0.22)  [Sauren  1993] 

Human  head:  15  mm  x  20  mm 
Skull  thickness:  10  mm 
BOUND  =  7.5  mm  +  2.0  A.air 
Xinc  =  BOUND  -  0.4  Xm 
cpinc  =  0°  or  45°  (see  Fig  14) 

Pi„c  =  1  Pa 
DPW  =  20 
ITS  =  T/20 


Fig  14:  (a,  top)  Incident  angle  cpinc  =  0°.  (b,  bottom)  <pinc  =  45°. 


Fig  15:  Four  locations  (A,  C,  F,  H)  along  inner  and  outer  skull 
surface  which  are  at  approximately  the  left  and  right  ear  locations 
near  each  side  of  the  skull. 


0ute*  Air  Aceustk  Pressure  v*.  Time  (nnmlm#  Surface  Ace-vitk  Pre«urev*.Time 


Fig  16:  Instantaneous  acoustic  pressure  waveforms  recorded  at  four  locations  (see  Fig  15). 
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Outer  Air  Acoustic  Intensity  vs.  Time 


Fig  17:  Instantaneous  acoustic  intensity  distributions  at  a  frequency  of  3  kHz  at  four  locations  (see  Fig  15). 
(top  panels)  (pinc  =  0”  and  (bottom  panels)  <pinc  =  45°. 


Fig  18:  The  Analyze  human  head  MRI  image  slice  is  segmented  to  derive  the 
outer  head  contour. 
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The  Analyze  human  head  model  contains  a  large  number  of  points  (order:  105).  Therefore,  this  model  was 
further  simplified  (order:  10)  while  maintaining  the  original  outer  geometry.  The  material  internal  to  the 
outer  3D  head  surface  was  constructed  of  skull  material  properties  to  demonstrate  that  it  was  feasible  to 
quantify  the  sound  pressure  amplitude  at  a  location  within  the  head  relative  to  the  air-borne  sound  pressure 
amplitude.  The  3D  FEA  model  is  shown  in  Fig  19. 


Parameters: 

Outer  medium:  air  (cair  =  340  m/s, 

X,air  =  0. 1 1 3  m,  pair  =  1 .2  kg/m3) 

Inner  medium:  water  (cwater =  1500  m/s, 

pwater=  1000  kg/m3) 

Human  skull  (p  =  1412  kg/m3, 

Young’s  modulus  =  6.5  GPa, 
Poisson’s  ratio  =  0.22) 

Human  diameter:  ~  0. 1 8  m 
BOUND  =  0.09  m  +  0.9  K*  =  0.192  m 
Xi„c  =  BOUND  -  0.4  K\x  =  -0.1467  m 
Pine  =  1  Pa 
ITS  =  TOO 


Fig  19:  The  3D  FEA  human  head  model. 

Admittedly  this  is  an  unrealistic  realization  of  the  head  but  one  that  can  be  evaluated  in  three  dimensions  and 
allows  for  the  evaluation  quantitatively  of  propagation  into  skull  material.  A  transient  3-kHz  one-cycle  sine 
wave  was  incident  from  air  onto  the  simple  human  head  model  and  acoustic  pressure  waveforms  were 
obtained  at  three  locations  (Fig  20).  The  instantaneous  acoustic  intensity  peak  (Eqn  8;  §4.5)  in  air  was 
determined  to  be  about  2.7  mW/m2  and  the  acoustic  intensity  in  the  skull  bone  was  about  1.24  pW/m2.  The 
acoustic  loss  (Eqn  9;  §4.5)  across  the  3D  skull  model  surface  was  estimated  to  be  approximately  32  dB,  quite 
consistent  with  theoretical  estimates  (33  dB). 


Fig  20:  Acoustic  pressure 
waveforms  at  three  selected 
locations  for  the  3D  FEA 
human  head  model. 


TIME 
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These  human  head  studies  identified  a  problem  that  we  knew  we  were  going  to  have:  how  to  represent 
graphically  and  quantitatively  a  3D  process?  This  work  with  the  simpler  Analyze  human  head  model 
provided  needed  insight  to  represent  a  3D  process.  Therefore,  while  the  work  continues  to  develop  a  3D 
human  head  from  the  more  complete  NIH  dataset,  a  procedure  was  initiated  to  represent  graphically  and 
quantitatively  a  3D  transient  process. 

5.3  Ray  Tracing 

A  ray  tracing  approach  to  represent  graphically  and  quantitatively  a  3D  transient  process  has  been  developed. 
Ray  tracing  is  a  common  procedure  by  which  wave  propagation  is  displayed.  For  example,  in  the  study  of 
optics,  rays  are  used  to  depict  the  path  or  paths  taken  as  a  light  wave  travels  through  a  lens.  However,  in 
optics,  the  eikonal  equation  can  be  solved  because  the  wavelength  is  assumed  to  be  zero  so  that  propagation 
laws  can  be  formulated  in  terms  of  geometry.  This  is  also  the  case  for  geometric  acoustics  and  is  often  used 
to  solve  acoustic  propagation  problems  in  the  ocean.  However,  for  the  case  of  an  acoustic  wave  that  has  a 
wavelength  comparable  to  the  object  onto  which  it  is  incident,  the  full  wave  equation  must  be  solved  because 
diffraction  needs  to  be  included  as  part  of  the  analysis.  Therefore,  ray  paths  need  to  be  deduced  from  the 
propagated  acoustic  wavefront.  The  transient  ANSYS  FEA  procedure  yields  data  from  which  the  propagated 
acoustic  wavefront  has  been  deduced.  From  the  propagated  acoustic  wavefront,  rays  have  been  calculated  by 
taking  the  normal  components  of  the  wavefront  surface  as  a  function  of  time.  The  transient  FEA  procedure 
yields  data  as  a  function  of  time  steps,  thus  allowing  the  wavefront  to  be  animated  as  a  function  of  time. 
Also,  from  these  data,  as  long  as  the  time  steps  are  sufficiently  small,  ray  paths  can  be  calculated.  Further, 
after  the  ray  paths  have  been  calculated,  the  acoustic  energy  can  be  determined  from  the  density  (number  of 
ray  paths  that  intersect  a  unit  volume)  of  ray  paths  in  a  particular  volume. 

Transient  (3-kHz  one-cycle  sine  wave)  ANSYS  FEA  procedures  were  performed  using  a  9-cm-radius 
hemisphere  of  known  propagation  speeds.  The  plane  wave  was  incident  in  air  (speed  =  340  m/s)  onto  the 
spherical  surface  of  the  hemisphere  that  had  propagation  speeds  consistent  with  water  (speed  =  1540  m/s). 
Wavefronts  were  deduced  from  a  time-domain  correlation  technique.  From  the  wavefronts,  ray  paths  were 

calculated  (Fig  21  depicts  ray 
paths).  The  ray  paths 
originated  in  air  and 
propagated  into  the 
hemisphere.  For  low-angle 
incidence,  the  ray  path 
directions  into  the  hemisphere 
were  consistent  with  those 
calculated  from  Snell’s  Law. 
However,  as  the  angle  of 
incidence  onto  the 
hemisphere  increased,  the  ray 
path  directions  into  the 
hemisphere  became 
progressively  greater  from 
those  expected  from  a  Snell’s 
Law  calculation.  We  are 
currently  investigating  the 
cause  of  this  phenomenon. 
We  suspect  that  the  ray  path 
directions  are  correct  because 
diffraction  phenomenon 
could  be  causing  this 
problem.  For  this  example, 
the  wavelength  in  air  is  1 1 .3 
cm  (radius  of  the  hemisphere  is  9  cm)  and  ka  is  5.0  (all  based  on  a  center  frequency  of  3  kHz  but  the 


3khz-hsphere3: 340/1 540  m/s  3khz-hsphere3: 340/1 540  m/s 
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Fig  21 :  Ray  paths  calculated  based  on  FEA  data. 
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frequency  content  of  the  one-cycle  wave  is  much  broader).  To  test  our  hypothesis,  we  conducted  the 
transient  ANSYS  FEA  analysis  at  10  kHz  using  the  9-cm-radius  hemisphere.  However,  at  the  high 
frequency  of  10  kHz,  we  hit  the  current  ANSYS  version  limits  on  elements/nodes  number  and  a  coarse 
meshing  is  used  in  the  FEA  model  that  compromised  the  computation  accuracy.  Comparing  the 
computational  results  and  the  theoretical  solutions  calculated  using  Snell’s  law,  up  to  0.7°  difference  is  found 
for  low-angle  incidence  while  up  to  10°  difference  is  found  for  larger-angle  incidence  (Fig  22).  The  next  step 
is  to  conduct  the  FEA  analysis  using  planar  surfaces  as  a  function  of  incident  angle.  A  well-documented  and 
verified  ray  tracing  procedure  will  provide  an  easy  to  understand  (and  visualize)  dynamic  process.  Also,  the 
budget  requests  funds  to  support  partially  the  unlimited  ANSYS  version. 


Figure  23  shows  the  wavefront  at  three  time  steps  as  the 
wavefront  propagates  through  the  hemisphere.  If  the 
shown  arrows  were  connected  at  each  wavefront  step,  the 
resultant  lines  would  depict  ray  paths. 

In  summary,  we  have  (1)  validated  finite-element  analysis 
(FEA)  general  processing  code  for  both  harmonic  and 
transient  solutions,  (2)  constructed,  simplified  and  verified 
transient  FEA  analyses  of  the  2D  NIH  human  head,  (3) 
demonstrated  FEA  analysis  of  the  3D  Analyze  human 
head,  and  (4)  developed  a  ray  tracing  approach  to 
graphically  and  quantitatively  represent  a  3D  transient 
process. 


Tima  Step  -  36 


Fig  23.  Three  time  steps  as  the  wavefront  propagates  through  the  hemisphere. 
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